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We present an efficient first-principies approach for caicuiating Fermi surface averages and spectrai 
properties of soiids, and use it to compute tfie low-fieid Hall coefficient of several cubic metals and 
the magnetic circular dichroism of iron. The first step is to perform a conventional first-principles 
calculation and store the low-lying Bloch functions evaluated on a uniform grid of fc-points in the 
Brillouin zone. We then map those states onto a set of maximally-localized Wannier functions, 
and evaluate the matrix elements of the Hamiltonian and the other needed operators between the 
Wannier orbitals, thus setting up an "exact tight-binding model." In this compact representation 
the fc-space quantities are evaluated inexpensively using a generalized Slater-Koster interpolation. 
Because of the strong localization of the Wannier orbitals in real space, the smoothness and accuracy 
of the fc-space interpolation increases rapidly with the number of grid points originally used to 
construct the Wannier functions. This allows fc-space integrals to be performed with ab-initio 
accuracy at low cost. In the Wannier representation, band gradients, effective masses, and other 
fc-derivatives needed for transport and optical coefficients can be evaluated analytically, producing 
numerically stable results even at band crossings and near weak avoided crossings. 



I. INTRODUCTION 

In electronic structure calculations for solids, the eval- 
uation of an observable requires integrating a periodic 
function in reciprocal spaceii We will distinguish between 
three kinds of properties: those where the integral is over 
the Brillouin zone (type-I), over the Fermi surface (type- 
II), and over an energy-difference isosurface (type-Ill). In 
many cases those integrals take, at T = 0, the form 
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Here is the Fermi level, are the eigenenergies of 
the one-electron states, and i^„,„(k) involves matrix el- 
ements of periodic operators which commute with the 
crystal translations. Ground-state properties such as 
the total energy, and dc response functions such as the 
Hall conductivity, are examples of the first and second 
type of property, respectively. The third type includes 
optical absorption in the dipolc approximation; other 
response and spectral functions can be expressed in a 
similar form.^ Eqs. II])-© are by no means exhaustive. 
While properties such as the elect ron-phonon interaction 
[described by matrix elements of the form Fnm (k, k 4- q) 
associated with phonon wavevector q] are not explicitly 
covered above, the methods discussed in this paper can 
be extended to handle such cases<^ 



In a practical calculation the continuous integral is re- 
placed by a summation over a finite number N of points 
in the Brillouin zone (BZ), 
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where Vcoii is the cell volume and i(j(k) are the fc-point 
weights that arise upon restricting the summation to the 
irreducible wedge of the BZ. For type-I properties of in- 
sulators, the integrand varies smoothly across the BZ 
and this summation converges rapidly with the number 
of sampled points.— In metals the BZ integral must be 
treated carefully, as the integrand is now discontinuous 
due to the partial filling of the energy bands. Prop- 
erties of type II and III pose the additional challenge 
of sampling isosurfaces in fc-space accurately and effi- 
ciently In all these other cases the rate of convergence 
is much slower than for ground state properties of in- 
sulators, and a very large number of fc'-points may be 
needed. This sampling problem severely limits the ef- 
ficiency and accuracy of ab-initio calculations for many 
properties. Examples of such difficulties abound even in 
the recent literature, and include the magnetocrystalline 
anisotropy of ferromagnets^ and optical absorption in hot 
liquid metalS)^ to name just a few. 

In this paper we describe a widely-applicable WF- 
bascd scheme for interpolating both the energy bands 
£„k and the matrix elements i^„,„(k). The method was 
used in Ref. to compute the anomalous Hall conductiv- 
ity of iron. Where possible, we have adopted a notation 
consistent with that work. Ref. @ dealt with a type-I 
problem (the quantity being integrated over the BZ was 
the Berry curvature of the occupied states), and here we 
extend the method to problems of type II and III. 

As an example of a type-II problem, we study the low- 
field classical Hall coefficient of several cubic metals. This 
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and other transport coefficients pose an additional chal- 
lenge to existing ab-initio methodologies: how to evaluate 
the first and possibly also the second fc-derivatives of the 
energy bands at the Fermi level. Early worki employed 
tight-binding (TB) parameterizations of ab-initio bands 
and the derivatives were calculated by numerical differ- 
entiation using the linear tetrahedron method. In other 
work, an analytic evaluation of the TB band gradients 
has been used to achieve improved numerical stability;^ 
but the second derivatives were still computed by finite 
differences. Other interpolation strategies, such as the 
SKW scheme^iiSiiii and spectral differentiatioufi^ have 
also been used. 

All previous interpolation schemes have one feature in 
common: the only information retained from the origi- 
nal ab-initio calculation is the set of energy eigenvalues 
on a grid of /s-points. Hence the information about the 
connectivity of the bands is lost, and the interpolation 
becomes unreliable or even unstable in the vicinity of 
band crossings, avoided crossings, and near-degeneracies. 
Moreover, retaining only the eigenenergies strongly re- 
stricts the type of matrix elements i^„,„(k), and hence 
observables, that can be computed. 

A more powerful interpolation scheme can be obtained 
by keeping one more piece of information, namely, the 
overlap matrices between the Bloch states at neighboring 
grid points as in Eq. Q below. It is perhaps not widely 
appreciated that the information about band connectiv- 
ity is encoded in those overlap matrices. Indeed, they 
are the key input for the WF-construction methodp^iii^ 
and the connectivity can be recovered from the Wannier 
representation of the bandstructure. Thus, not only do 
the Wannicr-interpolated bands reproduce the ab-initio 
bands with essentially no loss of accuracy, but their fc- 
derivatives can also be evaluated analytically. Like the 
SKW scheme, the present method is based on Fourier in- 
terpolation. Unlike SKWfiS however, it produces stable 
and reliable results even in the presence of band crossings 
and avoided crossings. 

Remarkably, a knowledge of the overlap matrices of 
Eq. ([5|) allows for the interpolation of properties that 
are not determined by the energy bands alone, but also 
depend on the position (or velocity) matrix elements. 
(More generally, any one-electron operator can be inter- 
polated if, in addition, its matrix elements between the 
WFs are tabulated.) As an example, we compute the 
magnetic circular dichroism of iron, a typc-III property. 

The paper is organized as follows. Section [III contains 
the methodological aspects of the work. We start by re- 
viewing the WF construction methods. We then describe 
the Wannier-intcrpolation strategy for a generic periodic 
operator. The interpolation of the velocity operator, as 
well as of band gradients and inverse effective masses, 
is discussed separately. We conclude Sec. HI] by present- 
ing an improved broadening scheme for performing the 
fc-space integrals. In Sec. IIIII we apply the technique to 
the low- field Hall effect of several cubic metals, and in 
Sec. IIVI to the magnetic circular dichroism of bcc Fe. In 



Sec.|V]wc provide a brief discussion and conclusion. The 
Appendix contains some convergence studies. 



II. WANNIER INTERPOLATION 
A. Construction of the Wannier functions 

Ab-initio calculations provide a certain number of low- 
lying Bloch eigenstates '(/'nq(r) = e''=' '"M„q(r) on a mesh 
of fc-points in the BZ, which we take to be uniform. We 
will denote points on this ''^ab-initio mesh" by q, to distin- 
guish them from arbitrary or interpolation points, which 
will be denoted by k. 

Consider a type-II (Fermi-surface) problem; two situa- 
tions may occur. The ffi'st one, which is seen in Pb for ex- 
ample, occurs when the Fermi level lies within an isolated 
group of AI bands, where by "isolated" we mean sepa- 
rated from all higher and lower bands by a gap through- 
out the BZ. In this case it is possible to construct a set 
of M WFs per unit cell spanning the Hilbert space of 
the isolated Bloch manifold. This can be done using the 
method of Marzari and VanderbiltH to obtain so-called 
maximally localized Wannier functions for that isolated 
group of bands. 

The second scenario occurs when the bands of inter- 
est are "entangled" with other bands. Then it is still 
possible, using the approach of Souza, Marzari, and 
Vanderbilt^ to construct a small number M of maxi- 
mally localized WFs which describe those bands exactly. 
The number M of WFs per cell is now to some extent 
an adjustable parameter. The first step is to identify the 
subspace of states of interest. Usually this is done by se- 
lecting the bands inside an energy window spanning from 
Emin to Eniax- (For typc-I and type-Ill problems E'min 
is normally in the gap below the lowest valence bands 
and the position of £^max depends on the problem, but 
is always above E{. For a type-II problem the only re- 
quirement is that -Eniin < E{ < i?inaxj ^-nd in practice the 
range is adjusted so that WFs with good localization and 
symmetry properties result.) The number of states 
within this window can vary from one q-point to another, 
and we require that M > Nq for all q, so that the space 
spanned by the WFs (the "projected space") can be cho- 
sen to contain as a subspace all the window states. In 
the method of Ref. [l3 a second (outer) energy window is 
used which encloses the previously defined (inner) win- 
dow. At each q, the M-dimensional projected space is a 
subspace of the A/'q-dimensional space of states contained 
in the outer window. For the special case of an isolated 
group of M bands it is natural to choose M = A'q = A/'q 

Only two pieces of information from the ab-initio cal- 
culation are needed as an input to the WF-gencration 
algorithm: the A/'q band-energy eigenvalues £„q, and the 
A/'q X A/'q+b overlap matrices between the cell-periodic 
Bloch eigenstates at neighboring points q and q + b, 

<5'„m(q, b) = (UnqlUm^q+b). (5) 
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The output consists of an Afq x M matrix U{q) for each 
q. These matrices relate the original set of Afq ab-initio 
Bloch eigenstates selected by the outer window to a new 
set of M orthonormal Bloch-like states 

771—1 

that vary smoothly with q. These states are labeled with 
a superscript (W) to indicate that the WFs arc obtained 
from them by a direct Fourier sum 

where the sum runs over a grid of Nq q-points and |nR) 
is the n-th Wannier function in the unit cell at R. 

Although the explicit construction of the WFs obvi- 
ously requires a knowledge of the |u„q)'s, only the eigen- 
values £nq and the overlaps S'(q, b) are needed to obtain 
the U{q) matrices. Retaining this minimal information 
from the ab-initio calculation is thus sufficient for many 
applications, including the ones presented in this work. 

An important object in what follows is the M x M 
Hamiltonian matrix in the projected subspace, 

= [uHq)H{q)U{q)]^^^, (8) 

where Hnmil) = ^nqSnm is a diagonal J\fq x A/'q matrix 
and i?(q) = e-^i'^TJe*'! ^ We diagonalize H^'^^q) by 
finding an M x M unitary matrix f/(q) such that 

C/t(q)if(W)(q)t/(q)^i7(H)(q)^ (g) 

where i?im(q) = SnqSnm- Thcn £',lq'' will be identical 
to the original ab-initio Snq for all bands inside the inner 
window. The corresponding Bloch states 

i<') = Ei""^^)t^™(q) (10) 

771 

will also coincide with the ab-initio eigenstates |M„q) in- 
side the inner window. We shall refer to a quantity with 
a (W) or (H) superscript as belonging to the Wannier or 
Hamiltonian gauge respectively. 

B. Wannier interpolation of a periodic operator 

The problem we pose for ourselves is the following one. 
Suppose we are given a periodic operator operator O, and 
we have computed at every q 

C7l777 (q) = (w77q|0(q)|Mr77q>, (11) 

its matrix elements between the A/'q ab-initio eigenstates 
in the outer energy window. How can we interpolate 



this matrix onto an arbitrary point k? We now show 
that this can be achieved once the matrices U (q) and the 
eigenvalues £nq (n — 1, . . . , A/'q) are known. Naturally, 
we can only expect the interpolation onto a given k to 
be meaningful for those matrix elements (n, to) for which 
both £nk and £mk fall within the inner window. 

We start by describing in Sec. IH B II the interpolation 
strategy as it applies to most ( "conventional" ) proper- 
ties. Transport and optical properties merit a separate 
discussion, given in Sec. Ill B 21 

1. Conventional properties 
By analogy with Eq. ^ , we define the M x M matrix 

= [Wt(q)(9(q)W(q)]_. (12) 
Next we find its Fourier sum 

Oiyj(R) = ^E^"^"'''^^™'(q)- (13) 

This operation is done once and for all for each of the 
A''o lattice vectors R lying in a supercell conjugate to the 
q-mesh. (If the sum is performed using a fast Fourier 
transform (FFT), the vectors R will be disposed in a 
parallelepipedal supercell.) Using Eq. ([7]) we recognize 
in oim' (R) the matrix element of O between WFs, 

OZH^) = {nO\d\mK). (14) 

In the above equations, the specification of the lattice 
vectors R can be left ambiguous with respect to supercell 
translations (R R + Rgup) since exp(iq- Rgup) = 1 for 
ah mesh points q, and thus O^^) (R + R-sup) = O'^'^^R). 
However, we now wish to perform the inverse (slow) 
Fourier transform 

a(r(k) = ^e^'^-^o(r(R), (15) 

R 

which yields the interpolation of Eq. ([T^ onto an arbi- 
trary point k. At this point the set of lattice vectors 
must be defined more precisely, since for points k not 
on the q-mesh exp(ik • Rgup) 7^ li and the smoothness 
of interpolation will depend on the choice of set. Us- 
ing the FFT parallelepipedal supercell, for example, is 
generally not optimal. Instead, one wants to choose lat- 
tice vectors lying inside the Wigner-Seitz (WS) super- 
cell centered on the origin^ii^ but the details may vary 
(e.g., sharing weights of R- vectors lying on the bound- 
ary of the WS supercell, or truncation to a sphere lying 
within the WS supercell). In practice the |Oim^(R)| de- 
cay exponentially with |R|, as expected if the WFs are 
exponentially localized, so the results should not be very 
sensitive to this choice, when using a sufficiently dense 
q-mesh. These remarks are illustrated in the Appendix. 
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FIG. 1: (Color online.) Wannier-interpolated bands of bcc Fe 
along F-H. The bands are color-coded according to the value 
of {Sz): red for spin up and blue for spin down. The energies 
are given in eV and the Fermi level is at eV. The vertical 
dashed lines indicate fc-points on the ab-mitio mesh used for 
constructing the WFs. 



The final step is to transform the matrix of Eq. pS]) 
from the Wannier to the Hamiltonian gauge. To find the 
required unitary matrix t/(k) we repeat the above steps 
ioT O = H to obtain H^'"'>{k.). The matrix C/(k) is then 
given by Eq. ([9]), with the replacement q — > k. Finally, 



OW(k) ^ J7t(k)0(W)(k)[/(k)^ 



(16) 



where M x M matrix products arc implied on the right- 
hand side. This solves the problem posed above at the 
beginning of Sec. IIIBI 

Once the WF matrix elements of both the operator of 
interest and the Hamiltonian are tabulated, the interpo- 
lation onto an arbitrary fc-point requires only inexpensive 
operations on small M x M matrices. When O = H, 
the present scheme reduces to Slater-Koster interpola- 
tion, with the maximally localized WFs playing the role 
of the TB basis orbitalsj^ 

Fig. 1 shows the interpolated band structure of bcc Fe 
along F-H, using the same WFs and computational de- 
tails as in Ref. [^. If one were to superimpose on this plot 
the energy bands obtained by performing an ab-initio cal- 
culation for a large number of points along the same line 
in fc-space, they would be essentially indistinguishable 
from the interpolated ones. Following Ref. [l^ we indi- 
cate with vertical dashed lines /c-points on the q-mesh 
used for constructing the WFs (an 8x8x8 grid in the 
full BZ). It is apparent that the Wannier interpolation 
procedure succeeds in resolving details on a scale much 
smaller than the spacing between those points. In par- 
ticular, the correct band connectivity is obtained. This 
means that spin-orbit-induced avoided crossings, for ex- 
ample, are never mistaken for actual crossings, no matter 
how weak the spin-orbit interaction. 

We note in passing that one could have formulated the 
problem at the beginning of Sec. Ill Bl somewhat differ- 
ently: rather than viewing the matrix elements of O and 



H between the Afq ab-initio Bloch states at each q as the 
basic ingredients of the method, we could have assigned 
that role to the matrix elements of the two operators be- 
tween the WFs. Even if the latter viewpoint is in some 
ways the more fundamental one, in practical implemen- 
tations one often obtains the Wannier matrix elements 
(ITil) via Eqs. ([ni)-(1131)- When doing so, the Wannier 
orbitals are never explicitly constructed. 



2. Transport and optical properties 

The treatment of transport and optical properties in 
crystals is more subtle. We will restrict our discussion to 
the electric-dipole approximation, where those properties 
depend on matrix elements of the velocity operator. The 
formulation of the previous subsection could in principle 
be used to interpolate the full velocity operator Uq, — 
~{i/h)[fa,H] {a = 1,2,3). Its matrix elements, as those 
of any other periodic operator, transform between the 
Wannier and Hamiltonian gauges according to Eq. pB]) 
(such a matrix will be called "gauge-covariant"— ) . They 
are given hy^ 



Unm,Q(k) = (V'nk|'Ca|'0mk) = 



Oka 



(17) 

However, when describing the dynamics of electrons in 
crystals it is convenient, from the points of view of both 
numerics and physics, to decompose the velocity operator 
into two terms.— This is achieved by taking da — d/dka 
outside the brackets on the right-hand side of Eq. pT|) 
and compensating for the extra terms that appear. After 
a few manipulations one obtains 



h dka 



4- (^?7ik ^nk)^Ti 

h 



where 



Ar. 



,(k) = i{Un-k\daU„ 



,.a(k) (18) 



(19) 



The first (second) term on the right-hand side of 
Eq. p8)) is diagonal (off-diagonal) in the band index. 
Note that neither is separately gauge-covariant. For ex- 
ample, the second one contains fc-derivatives acting on 
the eigenstates. According to Eq. pU]) . additional terms 
involving dU(\i)/dka will therefore appear when trans- 
forming between the Wannier and Hamiltonian gauges. 
Let us define, for every matrix object O, 



O 



(20) 



so that, by definition, O = O^^^ only for gauge- 
covariant objects. This notation will be used in the next 
section for expressing the intraband (diagonal) velocity 
matrix elements and the effective mass tensor, and in 
Sec. llVl for the intcrband (off-diagonal) velocity. 
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C. Band gradient and Hessian 

1. Notation 

We make use of the first and second fc-derivatives of 
the Hamiltonian matrix, 



Hn 



dka 



dkadkfs ' 



(21) 



(22) 



and define Hnm,a and H„m,af3 via Eq. (j20p as usual. We 
also define the first and second /c-derivatives of the band 
energy, 



Vnk,( 



1 „i< 



h dka 



1 g^gnk 



(23) 



(24) 



which have the interpretation of group velocity (ignoring 
Berry-curvature contributions) and inverse effective mass 
tensor, respectively. The strategy for interpolating these 
quantities is similar to the one developed in Ref. @ for 
the Berry curvature. We will again make extensive use 
of the antihermitian matrix 



T7(H) 

nrn^a 



defined in that work. 
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2. N on- degenerate hands 



First we consider the band-gradient velocity, Eq. 



In the Hamiltonian gauge H, 



(H) 



On On 



(H) 
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H, 

to /tQ. , 



/ivIq-'J, 



and hence 
Differentiating Eq. © with respect 



= S^™ + {i7mi?m+h.c.}, 



(26) 



where each object is an Af x M matrix and h.c. denotes 
the Hermitian conjugate. Because of the extra terms in 

curly brackets we have i/i^"* H^^^ and thus Ha, the 
first term on the right-hand side of Eq. p8)) . is not gauge- 
covariant. However, those extra terms do not contribute 
to the velocity; being the product of a diagonal matrix 
(^(H)^ with an antihermitian matrix {D^'a^), they only 



contain off-diagonal elements which cancel those in H ^ 
Thus 



;iv(") = = H 



T7(H) 

nn.a 



(27) 



Differentiating this equation yields the inverse effective 
mass tensor (l24t 



*2 (H) 



T7(H) 



U^H^^^U] +{u^H(J'^d0U + h.c.} 

J nn J nn 

{i?i«^<'+h.c.} . (28) 

1^ J nn 



(H) 

nn.a [3 



Unlike Eq. ([26|) . here the matrix in curly brackets has 
nonzero diagonal elements which contribute to /x^ 

Eqs. (P7|) - (P5)) are the desired expressions for the band 
derivatives, valid away from degeneracies and inside the 
inner energy window. They involve the M x M matrices 
U{k) calculated in Sec. [HBI] given by Eq. ([25|l. 

and H^^\ The last two involve fc-derivatives of 
Eq. (fT51) that can be taken analytically, i.e., 

H^ZU^) = E e*-^*i?.(nO|i/|mR> (29) 



R 



and 



hZIpO^) = E e"^-''{~RaRp){nO\H\mK). (30) 



R 



5. Discussion 

In order to interpret the above expressions it is illu- 
minating to introduce ||(/)n)), the n-th M-component col- 
umn vector of ||0n)) is an eigenvector of H^'^\ the 
Hamiltonian operator projected onto the WF space. We 
then recognize in Eq. (|27[) the Hellmann-Fcynman result 

^ivia^ = ((0„||i?Q^^||0n)), and in Eq. ((28|) the expression 
for the effective mass tensor in empirical TB theory.— 
Eq. (P5)) is the standard result from k • p perturbation 
theory, in terms of the TB states. It involves the oper- 
ator {l/h)H^\ which differs from the full velocity op- 
erator in that the position-operator-dependent terms are 
absent f^i^ We note that all the formulas given so far and 
in the rest of the paper remain valid when the ab-initio 
Hamiltonian contains non-local and spin-orbit terms. 

The advantage of this reformulation of k • p perturba- 
tion theory is that it is done strictly in terms of the small 
number M of M-dimensional states \\4>n}} at an arbitrary 
k, and yet it is exact within the inner energy window. 
In contrast, the formulation in terms of the original ab- 
initio states on the q-grid is considerably more expensive 
and usually entails a truncation crrorii^ 



4- Degeneracies 

While meaningful band derivatives can be defined via 
degenerate k • p perturbation theory even at degeneracy 
points in the BZji^ this is not possible when the only 
information available about the bandstructurc is a list of 
eigenenergies on a predetermined coarse fc-point grid. In 
that case, the information about the band connectivity 
is lost, and finite-difference estimates of the derivatives 
become ill-defined at points of degeneracy, which have to 
be carefully avoided.— 

In the present formulation the fc-gradients of the de- 
generate states are the eigenvalues of the submatrix 



H 



(H) 



fiiy 



(31) 



where the indices ^ and v run over the degenerate states 
only. We update the M x M matrix U by replacing the 
columns corresponding to those states with the rotated 
states that diagonalize The Hessian matrix can 



then be obtained from Eq. ([28|) using the updated U and 
a modified form of Eq. (P5)) . 




FIG. 2: Upper panel: Dispersion of the three p-like energy 
bands of Pb along the F-K direction, obtained by interpo- 
lating a non-relativistic ab-initio calculation. Lower panel; 
Inverse effective masses of those bands along the same direc- 
tion, calculated in two ways: from the interpolated eigenen- 
ergies on a regular mesh of points using a spline fit (circles), 
and from perturbation theory in the Wannier representation 
(solid lines). 



£)(H) 



T7(H) 
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o(H) _ ^(H) 
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(32) 



In cases of degeneracies (at band edges, for example) 
where some of the eigenvalues of the matrix are 
equal, a first-order treatment is inadequate and the cor- 
rect rotation between eigenstates needed to compute the 
Hessian is found by going to second order in degenerate 
perturbation theory. This amounts to diagonalizing the 
submatrix obtained from the right-hand side of Eq. (pS]) 
by replacing the subscripts nn therein by fiiy, and letting 
/i and v run over the the first-order-degenerate bands. 
The desired Hessian matrix elements are the eigenvalues 
of that submatrix. 

In our calculations we employ the first-order form of 
degenerate perturbation theory when two or more energy 
eigenvalues lie within 10~^ eV of each other and we sub- 
sequently use the second-order form if in addition two or 
more band gradients differ by less than 0.1 eV-A. 



with the PWSCF codo^ using density-functional theory 
in the local-density approximation, together with the 
plancwave-pseudopotential formalismFi Norm-conserving 
pseudopotentials are employed, and spin-orbit effects are 
included in Sec. llVl and Fig. [1] only. The WFs are gener- 
ated using the WANNIER90 code^ 

In Pb the lowest four valence bands crossing the Fermi 
level are separated everywhere in the Brillouin zone from 
higher bands, so that the original prescription of Marzari 
and Vanderbilti^ can be used to generate the correspond- 
ing maximally localized WFs. They are atom-centered 
and have sp^ character, with the axis of each orbital 
pointing towards a nearest neighbor. 

The inverse effective masses arc shown as solid lines 
in the lower panel of Fig. [21 For comparison, we also 
plot as circles in the lower panel the values obtained by 
fitting a spline function to the energy eigenvalues on a 
dense grid of fc-points along the path. We remark that 
whereas in the analytic method band crossings arc han- 
dled automatically, in order to obtain a smooth spline fit 
it was necessary to manually reorder the eigenvalues close 
to the band crossing. Once that is done, the agreement 
between the two is essentially perfect. 



5. Validation 

As an illustration, we used the above formulas to calcu- 
late the inverse effective mass of the three p-like valence 
bands of Pb along the F-K direction in fc-space. (We 
have chosen this example because it displays a threefold 
band-edge degeneracy at F and a band crossing between 
F and K.) 

In all the calculations in this work, lattice constants 
are taken from Ref. [2^ The Bloch states are obtained 



D. Adaptive broadening scheme for fc-space 
integration 

We conclude this section by discussing the evaluation 
of typc-I, II, and III integrals, Eqs. (dHS]). In order to ac- 
celerate the convergence of type-I integrals with respect 
to the number of sampling points in Eq. ^ , a broadening 
scheme can be usedj^i^ This amounts to replacing the 
step function in Eq. ([T|) with a Fcrmi-Dirac-likc smear- 



7 



ing function. In the case of type-I integrals, smearing is 
most important when relatively few sampling points are 
used, as tends to be the case whenever the evaluation of 
the integrand is expensivci^ If, however, the integrand 
is inexpensive, as is the case when using Wannier inter- 
polation, then it is possible to converge the BZ integral 
without resorting to smearing. For example, no smear- 
ing was used in Rcf . for integrating the Berry curvature 
over the occupied states of bcc Fe. 

Smearing plays a more fundamental role in integrals 
of types II and III: when replacing the BZ integral in 
Eqs. (OS]) by a grid summation, the (5-functions must be 
replaced by normalized functions with non-zero width, 
such as Gaussians. For example, in Eq. ([2]) one would 
replace ^(-E'f — £nk) by 



1 



2ttW 



cxp 



— {Ef — SnkY 

2W^ 



(33) 



Ideally the Gaussian width should be, for a given grid 
spacing Ak, comparable with the level spacing ASnk- 
The level spacing is however difficult to estimate, and the 
common practice is to set to a constant for all bands 
and fc-points. As a result, FS sheets arising from steep 
and flat bands are not described consistently. This is a 
serious disadvantage of broadening schemes with respect 
to the linear tetrahedron method, as discussed in Ref . [25l. 

This drawback is easy to remedy within the Wannier- 
interpolation method, since the band derivatives are 
readily available (Sec. IIICp . and can be used to esti- 
mate the level spacing. The simple estimate A£„k ^ 
l^^nk/SkjAfc suggests using a state-depending broaden- 
ing width 



Wn 



Ak 



(34) 



for type-I and type-II integrals (a is a dimensionless con- 
stant of the order of unity), and 



mk 



9k 



ak 



Afc 



(35) 



for type-Ill integrals. With this prescription W is no 
longer an independent adjustable parameter from Afc, 
guaranteeing that the Afc — > and W ^ limits are 
approached consistently. Several smearing functions be- 
yond a simple Gaussian have been propose d^^'^^ and can 
be used straight forwardly with the adaptive smearing 
scheme. For all of the calculations presented in this work 
we use the first-order Hcrmite polynomial scheme intro- 
duced by Methfessel and Paxtoni^ 

The above first-order adaptive smearing should be re- 
liable whenever the level-spacing is gradient-dominated. 
In practice we find that it works rather well even near 
critical points. This is illustrated in Fig.[3l where we show 
the density of states of diamond calculated using both the 
adaptive and conventional (fixed width) smearing, with 
a 50 X 50 X 50 interpolation mesh. When using a fixed 



Fixed Smearing (0.4eV) 
Adaptive Smearing 




Fixed Smearing (0.2eV) 
Adaptive Smearing 



12 14 
Energy (eV) 



FIG. 3: (Color online.) Density of states of bulk diamond 
calculated in tlie range 8-18eV using the conventional Gaus- 
sian broadening approach (light thin lines) with a fixed width 
of 0.4 eV in the upper panel and 0.2 eV in the lower panel, 
versus the adaptive broadening approach (dark thick lines). 
The inset shows the density of states in the full valence band 
range computed with the adaptive broadening approach. 



width of = 0.4 eV, the sharp van-Hove features are 
not well described. Reducing it to 0.2 eV improves the 
situation for some of them, but introduces spurious os- 
cillations whenever the level-spacing becomes larger than 
W. With the adaptive scheme such oscillations do not 
occur for a sensible choice of a, and the sharp features 
are well-described. We have used a = 1.0, but find the 
results to be quite robust for 0.8 < a < 1.3. 



III. LOW-FIELD HALL COEFFICIENT OF 
CUBIC METALS 

As a first benchmark application we compute the "clas- 
sical" low-field Hall coefficient of cubic metals, which is 
given by 



and 



dk 



(27r) 



dk 

(27r)3 



3 "'"nkV^k,: 



5/ 
d£ 



(36) 



(37) 



'rik 



a/ 

d£ 



(■^nk.xMnk,; 



(38) 



(For the systems studied in this Section, which are non- 
ferromagnetic and do not include the spin-orbit inter- 
action, the sum over spin-degenerate bands will give 
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TABLE I: Hall coefficient Rh, in units of 10""m^C"\ Ref- 
erences to the experimental data can be found in Ref. [27l . 





This work 


Ref. 2 


Ref. JO 


Experiment 


Li 


-12.7 


-12.8 


-12.4 


-15.5 


Al 


-2.5 


-1.7 


-3.4 


-3.43 


Cu 


-4.9 


-5.2 


-4.9 


-5.17 


Pd 


-11.9 


-6.0 


-17 


-7.6 



rise to factors of two, which are not written explicitly.) 
a^x is the longitudinal conductivity, and and (Jxy.z = 
do'xy/dBz, where axy is the antisymmetric (Hall) con- 
ductivity. The above expressions are obtained from a 
Bloch-Boltzmann description of transport; for a detailed 
derivation, see Ref. [l^- We note that Eq. ((38|) is written 
in a form which is specific to cubic metals. The quantities 
Vnk.a and /i„k,Q/3 are given by Eqs. (|23p and ((24l) . /(£) is 
the Fermi-Dirac distribution function, and < is the 
electron charge. At low temperatures {—df/d£) tends to 
5{£ — E{), and both axx and (Jxy.z become FS integrals 
of the form of Eq. ([2]). In the constant relaxation-time 
approximation r,ik drops out from Eq. (|36p so that Rn is 
completely specified by the first and second band deriva- 
tives at Ef. 

Calculations were done for Li, Al, Cu and Pd. Unlike 
Pb, in these metals the set of bands crossing the Fermi 
level is not isolated. Therefore, in order to generate max- 
imally localized WFs we first used the disentanglement 
procedure summarized in Sec. Ill Al to obtain an optimal 
projected space. The number of bands contained therein 
must be at least equal to the number of bands crossing 
the Fermi level. However, it is often desirable to extract a 
somewhat larger projected space if this produces a more 
symmetric set of Wannier functions. 

For lithium we obtained four atom-centered WFs per 
primitive cell with sp'^ character. For aluminum we ex- 
tracted a nine-dimensional projected subspace. The re- 
sulting WFs are atom-centered, but have no distinct sym- 
metry characteristics. For Cu and Pd we used seven 
WFs: five with d character on atom centers, and two 
with s character located at the tetrahedral interstitial 
sites. These have been previously described for Cu in 
Ref. B 

The computational details are the same as in Section 
III C 51 We obtain the self-consistent ground state us- 
ing a 16 X 16 X 16 Monkhorst-Pack mesh of fc-points and a 
fictitious Fermi smearing^! of 0.02 Ry for the Brillouin- 
zonc integration. We use the local density approxima- 
tion; for the materials studied we find the results are 
not significantly altered by using a GGA such as PBEi^ 
To compute the Hall coefficient we use an ab-initio grid 
of 12 X 12 X 12 q-points and obtain the required quantities 
on a uniform mesh of 60 x 60 x 60 k-points. We implement 
an adaptive mesh refinement scheme in which we iden- 
tify those points of the /c-space mesh at which at least 
one band lies with 0.5eV of the Fermi energy and obtain 
the required quantities on a 7 x 7 x 7 submesh spanning 



the original cell associated with this mesh point. We 
find these parameters give converged values of the Hall 
coefficient for the four metals studied. This is particu- 
larly reassuring in the case of Pd, where previous tech- 
niques encountered difficulties because of the occurrence 
of bands crossings near the Fermi surface 

The results for the Hall coefficient i?H are compiled in 
Table [D For Li, Al, and Cu wc find excellent agreement 
with previous calculations based on empirical TB fitting 
to ab-initio bands,— as well as ab-initio calculation com- 
bined with SKW interpolation^^ The case of Pd is more 
delicate as i?H depends critically on the details of the ab- 
initio calculation. For example a shift upwards (down- 
wards) in the Fermi level of 0.002Ry causes a change of 
—3 (+2) xlO~^^m'^C~^. In view of this wc refine the po- 
sition of the Fermi level in a final non-self-consistcnt step 
by interpolating the band energies and gradients onto a 
60 X 60 X 60 k-mesh and using the adaptive broadening 
scheme. Our converged value of Ru is intermediate be- 
tween the two previously computed values, and shows a 
relatively large discrepancy with experiment. Previous 
authors have suggestedii^ that it maybe necessary to go 
beyond the constant relaxation time approximation to 
give a good description of the transport properties of Pd. 

IV. MAGNETIC CIRCULAR DICHROISM 

Magneto-optical effects in ferromagnets result from a 
combination of exchange splitting and spin-orbit coupling 
(SOC),^i^ The former breaks time-reversal (TR) in the 
spin channel, and the latter transmits the TR-breaking 
to the orbital motion of the electrons, endowing the op- 
tical conductivity tensor with an antisymmetric compo- 
nent. The simplest such effect to evaluate is magnetic 
circular dichroism (MCD), the difference in absorption 
between left- and right-circularly-polarizcd light, and we 
have chosen it for illustrative purposes. It is given by 
the imaginary part of the antisymmetric conductivity, 

A. Evaluation of the Kubo formula 

Ab-initio calculations of magneto-optical effects de- 
mand high accuracy and dense /c-space sampling. The 
spin-orbit interaction is typically a small perturbation 
on top of the much larger exchange splitting, and the 
modifications that it produces on the electronic structure 
(both in the energy bands and in the matrix elements) 
are subtly and strongly dependent on fc-point and band 
index. 

(2) 

The conductivity ^^(w) is evaluated from the Kubo 
formula of linear-response theory in the electric-dipole 
approximation.— The needed ingredients are the energy 
eigenvalues of the states involved in the optical transi- 
tions and the transition matrix elements. We will eval- 
uate the interband contribution to the magneto-optical 
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absorption using Eq. ([18]) for the electric-dipolc transi- 
tion matrix elements, where it is now understood that 
all Bloch functions |Mnk) are spinors determined from 
a Hamiltonian that includes the spin-orbit interaction. 
One finds 

9 occ empty ,, 

n rn 

X [S{UJ - Ujnn) - S{UJ + Umn)] , (39) 

where hUmn = £m - £n- 

Eq. ([M)) is a type-Ill integral of the form of Eq. ([3]). 
When evaluating it by Wannier interpolation it must be 
kept in mind that the Wannier-derived bands reproduce 
the ab-initio ones only inside the inner energy window, 
and therefore its range must be adjusted according to 
the maximum desired absorption frequency. The matrix 
elements Anm.a arc to be evaluated in the Hamiltonian 

gauge, and the interpolation of Anm,a is based on the 
two relations^ 

4"^=X"^+*^i"^ (40) 

and 

4r.L(k) = E^''''''(0"l''"|I^'")' (41) 

R 

where D^^^ is given by Eq. jig) and A^f ' and ^L^^ are 
related by Eq. ^U^. Inserting Eq. (gOl) into Im(. . .) in 
Eq. dsn), we find 

^ Im [dI^I^D^^Ip) . (42) 

(2) 

The contributions to o'^ q^(w) from the three terms on 
the right-hand side will be denoted as A-A, D-A, and 
D-D, respectively. 

B. Results for bcc Fe 

Unlike the calculations presented earlier in the paper, 
to calculate the MCD spectrum we have used relativis- 
tic pseudopotentials which explicitly include spin-orbit 
effects.— The computational details for the ab-initio cal- 
culation, WF-generation, and treatment of the spin-orbit 
interaction are the same as in Ref. 0. The integral in 
Eq. (|39|) was evaluated on a uniform k-mesh containing 
125 X 125 X 125 points using the adaptive broadening 
scheme, which we find to be essential for resolving the 
fine details in the MCD spectrum. The spin magneti- 

(2) 

zation is along the z-axis, so that o'^ ^j^('-l') is the only 
independent non-zero component. 




Energy (eV) 

FIG. 4: Magnetic circular dichroism spectrum of bcc iron. 
The calculated spectrum (solid lines) is compared with the 
experimental spectrum from Ref. [s^ as reproduced in Ref. [S^ 
(open circles). 




Energy (eV) 

FIG. 5: (Color online.) Upper panel: decomposition of the 
magnetic circular dichroism spectrum into the three terms de- 
fined by the Wannier-interpolation procedure. Lower panel: 
cumulative anomalous Hall conductivity (AHC) versus lj, de- 
fined as the contribution to the AHC from frequencies higher 
than uj in Eq. (gSJ . 



It is conventional to plot the MCD spectrum as 

(2) rm 

ujap^\^y{uj), and adopt this convention in Fig. 01 Our 
calculated spectrum for bcc Fe is in excellent agree- 
ment with the one computed in Ref. [s^ using a com- 
pletely different electronic structure method. (Previ- 
ous calculations of magneto-optical effects have most 
commonly used all-electron methods.) As remarked in 
Ref. HO, this level of agreement between two different 
calculations is non-trivial when it comes to the spin- 

(2) 

orbit-induced ox j.y{^)- It provides a strong validation 
of the Wannier-interpolation scheme combined with the 
pseudopotential-planewave method as a viable way of 
computing magneto-optical effects. 

The upper panel of Fig. O shows the decomposition of 
the calculated MCD spectrum into the three contribu- 
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tions {A-A, D-A, and D-D) defined by the Wannier- 
intcrpolation procedure, as discussed following Eq. (j42|l . 
It is clear that the D-D contribution tends to dominate 
the spectrum in the frequency range from to 7 eV, es- 
pecially at the lowest frequencies. For frequencies above 
7 eV (not shown), the A-A and D-A terms become sig- 
nificant. 

The interband MCD spectrum (j)^]j.y{uj) is related 
to the Karplus-Luttingcr anomalous Hall conductivity^ 
(AHC) o'ALy(O) by the Kramers-Kronig relation 

O />OC -1 

^A,U(0) = -/ -/Zyi^)'^- (43) 

In the lower panel of Fig. [5] we show the cumulative AHC 
versus uj, defined as the contribution to the AHC from 
frequencies higher than lo in Eq. (|43p . In practice we use 
as the upper frequency limit in Eq. (j43|) the difference 
from the Fermi energy to the top of the inner energy 
window (18 eV). It is clear that the AHC is completely 
dominated by the low-frequency contributions below ^ 
5.5 eV. 

It can be shown that applying the transformation (I43p 
separately to the D-D term of the MCD spectrum yields 
the D-D term of the AHC, as defined in Ref. i. This 
explains the intriguing result that more than 99% of the 
anomalous Hall conductivity can be recovered from the 
D-D term alone.— This is a consequence of (i) the low- 
frequency part of the spectrum being weighted more in 
the integral as a result of the l/w factor in the integrand, 
and (ii) the D-D term overwhelming the other two at 
very low frequencies. 

V. CONCLUSIONS 

We have presented a Wannier-interpolation scheme 
to compute efficiently and accurately Fermi-surface and 
spectral properties from first principles. As an example 
of the former we computed the low-field Hall conductiv- 
ity for several cubic metals. As an example of the latter 
we calculated the magnetic circular dichroism spectrum 
of bcc Fe. 

The scheme naturally resolves a number of difficul- 
ties which have plagued existing interpolation schemes. 
Firstly, by preserving the information about band con- 
nectivity, band crossings and avoided crossings are 
treated correctly. In addition, the evaluation of the ve- 
locity matrix elements needed to compute both the Hall 
coefficient and the MCD spectrum can be done analyt- 
ically in the Wannier representation. Furthermore, the 
scheme does not become any more complex upon inclu- 
sion of the spin-orbit interaction in the Hamiltonian. In 
particular, there are no additional contributions to the 
velocity matrix elements; all the spin-orbit-related cor- 
rections are contained in the spinor WFs. Also, the 
Wannier-interpolation scheme is decoupled from the par- 
ticular choice of basis set used for performing the original 



ab-initio calculation, nor does it depend on the specific 
level of single-particle theory. As such, the calculation of 
a given property can be implemented in a universal way 
inside the Wannier module, which can then be interfaced 
with any desired electronic structure code. 

The appeal of the present approach is that it combines 
the simplicity of a tight-binding-like scheme with the 
power and accuracy of ab-initio methods. Most impor- 
tantly, it allows operators other than the Hamiltonian to 
be interpolated in the same manner as the Slater-Koster 
interpolation of energy bands. As such, it can be ap- 
plied to a wide variety of problems in condensed mat- 
ter physics. It should be particularly useful for study- 
ing metallic systems. A number of properties of met- 
als remain extremely challenging to compute from first- 
principles, as a result of difficulties in sampling the Fermi 
surface with sufficient accuracy. Wannier interpolation 
provides an elegant and powerful framework for investi- 
gating such problems with ab-initio techniques. 
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APPENDIX: CONVERGENCE PROPERTIES OF 
THE INTERPOLATION SCHEME 

For a given operator O the agreement inside the inner 
energy window between Onm (k) obtained by Wannier 
interpolation and OnmO^) calculated using a full first- 
principles calculation is determined by iVg, the number of 
points in the q-grid. The resulting WFs are periodic over 
the the conjugate real-space supercell spanning A^o unit 
cells. For any finite Nq there is some overlap between a 
WF and its neighboring periodic images, which affects 
the matrix ©(R). It is generally accepted that WFs 
decay exponentially; numerical studies have confirmed 
this for several materials;^ and recently there has been 
a claim of a formal proof for multiband time-reversal- 
invariant insulatorsj^ The error in ©(R), and therefore 
in the interpolation, should accordingly also decrease ex- 
ponentially beyond some supercell size. 

We report numerical tests for two cases: the isolated 
set of four valence bands in Pb, and the low-lying bands 
of Li, using the same WFs as in Sees. Ill C 51 and IIIIl re- 
spectively. The band energies are computed via both 
Wannier interpolation and non-self-consistcnt diagonal- 
ization of the planewave Hamiltonian on a 200 x 200 x 200 
BZ grid. For Li we collect data from the bottom of the 
inner energy window to 0.5eV below the top of the inner 
energy window; points close to the top of the inner win- 
dow may show larger discrepancies, as they result from 
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Linear dimension of the ab-initio grid 



FIG. 6: Convergence of the Wannier interpolated band ener- 

1 /3 

gies as a function of the Unear dimensions A'^p of the ab-mitio 
q-point grid. We plot the maximum error (squares) and mean 
absolute error (circles), where the error is the difference be- 
tween the Wannier interpolated band energy and the value 
obtained from a full non-self consistent diagonalization of the 
planewave Hamiltonian. The lines are linear fits to the points 
with A^o''^ > 8. 

an interpolation between q-points inside and outside the 
inner window. Fig. [S] shows several measures of the dif- 
ference in the energies as a function of N^^^. In both 



cases we find that the error decreases exponentially for 
-^0 ^10- It is particularly reassuring that this occurs in 
Li, since the decay properties of disentangled WFs has 
yet to be investigated thoroughly, and they probably fall 
outside the scope of existing formal proofs of exponential 
decay. 

Finally, we examine the optimal choice of supercell in 
which to define the set of lattice vectors R for the Fourier 
transform in Eq. (jlSp . To illustrate the discussion in- 
troduced earlier in the vicinity of Eq. (jisp . we compare 
the results for parallclcpipcdal and Wigner-Seitz super- 
cells. In the upper part of Fig. [7] we compare the inter- 
polated energy bands for a 4 x 4 x 4 grid of q-points. For 
such a sparse q-grid the interpolated bands do not agree 
precisely with the exact ab-initio bands from a non-self- 
consistent diagonalization of the planewave Hamiltonian; 
this is most noticeable in the deviation of the curvature 
of the three upper bands between K and T. However, it is 
clear that the Wigncr-Seitz supercell yields significantly 
better results than the parallelepipedal cell. This is most 
clear for the upper band from L to F, which displays large 
oscillations for the parallelepipedal cell. In the lower por- 
tion of Fig. [3 we show several measures of the error in the 

1/3 

interpolated bands as a function of A^q ; for any given q- 
grid the Wigner-Seitz cell gives the more accurate results. 
The superiority of the Wigner-Seitz choice can be easily 
understood, as it ensures the largest minimum distance 
between a WF and its periodic images. 
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FIG. 7: (Color online.) Comparison of interpolated band- 
energies for Pb obtained using sets of lattice vectors defined 
within Wigner-Seitz (WS) and parallelepipedal (P) supercells. 
Upper figure: Energy bands interpolated using a 4 x 4 x 4 
q-point grid (WS cell - solid lines, P cell - dashed lines). 
Full ab-initio results are shown as open circles. Lower figure; 
Convergence of the Wannier-interpolated band energies for 
the two supercells as a function of the linear dimensions A^q''"^ 
of the q-point grid. Details as in Fig. [6] 



